% reads edi files 
% rotate the tensors to a fixed angle
% projects the station locations to a profile
% returns the res,phase & errors of a single mode TE /TM
% latest date 1.8.4
% problmes - works only with procmt style edi files



function[Prof] = resprof(),

ProcDef.equations=[4 1 2;5 2 1];
ProcDef.neq=2;
ProcDef.rreqn = [4 1 2 6 7;5 2 1 7 6];
rotate = 0;
comp = 1; % 1 xy and 2 for yx    %<------------- Watch out for correct component !
%path1 = 'C:\Manoj\phd\finalxy\'; %<------------- Watch out for correct component !
path1 = 'C:\Manoj\PhD\PROCMT\';
%path1 = 'C:\Manoj\ngri-dell01\data\spectra\kutch\sn\';

DataCode = ['FILE:' path1 '_COMP:_' num2str(comp)];

minfrq = 1/200;
maxfrq = 200;

nfiles = 1;
%path1 = 'C:\Manoj\ngri-dell01\data\work\';
%path1 = 'C:\Manoj\ngri-dell01\data\spectra\sgt2\';

%path1 = 'E:\Manoj\data\spectra\sgt2\';
%fid = fopen('c:\manoj\phd\list.txt','rt'); %<----------------------------
dirs = dir(path1);
%S = fgetl(fid);  %<----------------------------
%while  S ~= -1,  %<----------------------------
%    S(S == ' ') = [];  %<----------------------------
%    S(end) = []; %<----------------------------
  for i = 3:length(dirs),
      fp = fopen([path1 dirs(i).name],'rt');
    if fp == -1,
        disp('Error in opening a file');
    %    disp([S S]); %<----------------------------
        disp([path1 dirs(i).name]);
    else,
    nfiles = nfiles +1;   
    [SPM,er]=edi_in(fp);
    ProcDef.nfrq=length(SPM.spectra);
    ProcDef.TLFreq1=SPM.head.frq; 
    ProcDef.file1=dirs(i).name;
    nchan = length(SPM.spectra(1).data);
    tf = TransFn(SPM,ProcDef);
    if rotate ~= 360 & rotate ~= 0,
        tf = rottf(tf,rotate); % rotate the tf
    end;
    [r,p,er,ep]=compresph(tf,ProcDef);
    
    [nfrq,tt] = sort(SPM.head.frq);
    Prof(nfiles-1).res = r(comp,tt);
    Prof(nfiles-1).reser = er(comp,tt);
    Prof(nfiles-1).ph = p(comp,tt);
    Prof(nfiles-1).pher = ep(comp,tt);
    Prof(nfiles-1).file = dirs(i).name;
    Prof(nfiles-1).site = SPM.head.site;
    Prof(nfiles-1).frq = nfrq;
    Prof(nfiles-1).Lat = SPM.head.Lat;
    Prof(nfiles-1).Log = SPM.head.Log;
    disp(dirs(i).name);
    %disp(ProcDef.nfrq);
    %semilogy(ones([1,ProcDef.nfrq])*nfiles,ProcDef.TLFreq1,'r.');
%    semilogy(ones([1,ProcDef.nfrq]),ProcDef.TLFreq1,'r.');
    %hold on;
    x(nfiles-1) = SPM.head.Log;
    y(nfiles-1) = SPM.head.Lat;
end;
    %S = fgetl(fid);%<------------------------------
    
end;

%fclose(fid);   %    <-------------
figure(2);

% for projection of sites to a line ->
[p,s] = polyfit(x,y,1);  % poly fit degree 1 (linear regress)
ymodel = polyval([p(1) 0],x);
a = [x(end);ymodel(end)]; % ponit through which the fitted line pass
% this should also be possible to get the desired line through GUI
P = (a*a.')/(a.'*a); %tranform matrix
U = P*[x;y-p(2)];  %check this when y intercept > 0
Prj_x = U(1,:); % projected X
Prj_y = U(2,:)+p(2); % projected Y

plot(x,y,'b*');
hold on;
plot(Prj_x,Prj_y,'r.');

for i=1:nfiles-1,
    Prof(i).Proj_Log = Prj_x(i);
    Prof(i).Proj_Lat = Prj_y(i);
end;

Prof = frequency_check(Prof,minfrq,maxfrq)
Prof = location_check(Prof);
Prof(1).DataCode = DataCode;
%end; % <-------------------------------


function[tf] = rottf(tf,rotate)

ang = rotate*pi/180;
b = [cos(ang),sin(ang);-sin(ang),cos(ang)];
for i = 1:length(tf),
    a = [tf(i,3),tf(i,1);tf(i,2),tf(i,4)];
    c = a*b;
    tf(i,1) = c(1,2);
    tf(i,2) = c(2,1);
    tf(i,3) = c(1,1);
    tf(i,4) = c(2,2);
end;


 function[Prof] = frequency_check(Prof,minfrq,maxfrq)

     for i = 1:length(Prof),
        L = Prof(i).frq > minfrq & Prof(i).frq < maxfrq;     
        Prof(i).frq(~L) = [];
        Prof(i).res(~L) = [];
        Prof(i).reser(~L) = [];
        Prof(i).ph(~L) = [];
        Prof(i).pher(~L) = [];
     end;

master = Prof(1).frq;
for i = 2:length(Prof),
for j=1:length(Prof(i).frq), 
    L = fix(log10(master)*100) == fix(log10(Prof(i).frq(j))*100); %To ingnor minor frequencies mismatch 
 % L = master==Prof(i).frq(j);
    if sum(L) == 0,
        master(end+1) = Prof(i).frq(j);
    end;
end;
end; 
master = sort(master);
for i=1:length(Prof),
    frq = Prof(i).frq;
       
    for j = 1:length(Prof(i).frq),
       array(j) = find(fix(log10(master)*100) == fix(log10(Prof(i).frq(j))*100));
 %   array(j) = find(master==Prof(i).frq(j));
    end;
    L = Prof(i).frq > minfrq & Prof(i).frq < maxfrq;     
    Prof(i).arr_frq = ones([1,length(master)])*NaN;
    Prof(i).arr_frq(array) = Prof(i).frq;
    Prof(i).arr_res = ones([1,length(master)])*NaN;
    Prof(i).arr_res(array) = Prof(i).res(L);
    Prof(i).arr_reser = ones([1,length(master)])*NaN;
    Prof(i).arr_reser(array) = Prof(i).reser(L);
    Prof(i).arr_ph = ones([1,length(master)])*NaN;
    Prof(i).arr_ph(array) = Prof(i).ph(L);
    Prof(i).arr_pher = ones([1,length(master)])*NaN;
    Prof(i).arr_pher(array) = Prof(i).pher(L);
    Prof(i).master_frq = master;
    clear array;
    
end;



function[Prof] = location_check(Prof)

for i = 1:length(Prof),
    x(i) = Prof(i).Proj_Log;
    y(i)= Prof(i).Proj_Lat;
end;

magn = sqrt(x.^2 + y.^2);  % this works in latitude ranges of India !!!
[magn,it] = sort(magn);
it = flipud(it'); % jut to keep it in descending order (there may be a problem if it is in | direction (not -  dir)

for i = 1:length(Prof), % there may me an elegant way of doing it. but i just do not know !
    NProf(i).res = Prof(it(i)).res;
    NProf(i).reser = Prof(it(i)).reser;
    NProf(i).ph = Prof(it(i)).ph;
    NProf(i).pher = Prof(it(i)).pher;
	 NProf(i).file = Prof(it(i)).file;
    NProf(i).site = Prof(it(i)).site;
    NProf(i).frq = Prof(it(i)).frq;
    NProf(i).Lat = Prof(it(i)).Lat;
    NProf(i).Log = Prof(it(i)).Log;
    NProf(i).Proj_Log = Prof(it(i)).Proj_Log;
    NProf(i).Proj_Lat = Prof(it(i)).Proj_Lat;
    NProf(i).arr_frq = Prof(it(i)).arr_frq;
    NProf(i).arr_res = Prof(it(i)).arr_res;
    NProf(i).arr_reser = Prof(it(i)).arr_reser;
    NProf(i).arr_ph = Prof(it(i)).arr_ph;
    NProf(i).arr_pher = Prof(it(i)).arr_pher;
    NProf(i).master = Prof(it(i)).master_frq;
    fprintf('%d, %d, %d \n',it(i),i,length(NProf));
end;


Prof = NProf;

Prof(1).dist=0;
for i = 2:length(Prof),
   Prof(i).dist = lat2km(Prof(1).Lat,Prof(i).Lat,Prof(1).Log,Prof(i).Log);
	%Prof(i).dist = deg2km(dist([Prof(1).Lat,Prof(1).Log],[Prof(i).Lat,Prof(i).Log]'));
   
end;
